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ROTATIONAL  DEFORMATION  OF  THE  EARTH 
AND  MAJOR  PLANETS 


INTRODUCTION 

The  purpose  of  this  paper  is  to  ascertain  the  deformations  sustained  by  a rotating 
self-gravitating  massive  body  when  hydrostatic  and  elastic  forces  are  present  within  its 
interior.  There  are  many  advantages  to  be  gained  by  studies  of  this  type  of  deformation; 
for  example:  (1)  better  determination  of  the  shape  of  the  geoid  in  conjunction  with 
accurate  gravity  measurements,  (2)  better  knowledge  of  the  Earth’s  bodily  tides  and  of 
the  Love  numbers,  (3)  effects  of  crustal  loading  on  the  global  ocean  tides,  and  (4)  deter- 
mination of  the  gravity  field  of  a planet. 

The  Earth  and/or  any  other  planet  are  considered  to  be  initially  in  hydrostatic  equil- 
ibrium under  the  influence  of  self-gravitatation  and  of  their  rotational  motion.  This  initial 
state  of  deformation  will  be  taken  as  the  reference  state. 

The  following  relationships  will  be  valid  between  the  stress  field  Tt]  and  the  hydrostatic 
pressure  p0: 


_ . dPo 

Tij  fa  ~ Po#0,i 


(1) 


Here  the  6's  are  the  Kronecker  deltas,  and  pQ  is  the  mass  density  of  the  reference  state. 
g0  i are  the  components  of  gravity  for  the  same  state  and  are  obtainable  from  the  gravity 
potential  as 


Bo.i 


avo 

dx. 


(2) 


This  potential  is  the  sum  of  two  terms:  the  gravitational  potential  Vj  and  the  rotational 
potential 


Vc  =“  ou2r2sin20. 


(3) 


* 


where  co  is  the  rotational  velocity,  r is  the  radius  vector  measured  from  the  center  of  gravity, 
and  0 is  the  colatitude.  V0  satisfies  the  Poisson  equation 


V2V0  = -4irGp0  + 2 to2, 
where  G is  the  gravitational  constant. 
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In  previous  publications  - Lanzano  1 1 1 and  Kopal  and  Lanzano  1 21  - we  used 
equation  (1)  to  ascertain  how  within  a massive  rotating  configuration  the  equipotential  sur- 
faces which  were  originally  of  spherical  shape  are  deformed  into  spheroids  having  an  equation 
of  the  sort 


— = 1 + 
a 


oo 

X (~)P2J  (cot°y 

i- o 


(5) 


Here  the  parameter  a denotes  the  mean  radius  of  the  spheroid,  a,  corresponds  to  the  outer- 
most equipotential  surface,  and  the  P's.  are,  as  usual,  the  Legendre  polynomials.  The  defor- 
mation coefficients  were  in  turn  expressed  as  power  series  of  the  dimensionless  parameter 


Q = 


co2af 

3 Cm, 


(6) 


where  m,  is  the  mass  of  the  total  configuration.  This  parameter  essentially  represents  the 
ratio  between  centrifugal  and  gravitational  accelerations.  The  accepted  value  of  q for  the 
Earth  is  0.00115,  for  Jupiter  0.02800-1,  and  for  Saturn  0.0-17207.  In  the  aforementioned 
publications,  we  developed  an  equilibrium  theory  valid  up  to  the  third  power  of  q;  we 
limited  ourselves  to  the  following  expansions: 


fo  - 

+ <13f03 

f 2 = ‘if 21 

+ <jV22 

+ Q3f 23 

u = 

*l'*f  42 

+ 1*f43 

f 6 = 

f'V6  3 

(7) 

The  unknown  functions  f^2  and  were  eliminated  via  mass  conservation  considerations 
and  found  to  be  related  to  f2 x and  f22  as  follows: 


fo2  "5^2?*  /"os  ~ ~ 6 ^21  ^22  + 21  ^2?)  • 


(«> 


The  other  six  unknown  deformation  coefficients  (i.e.,  f2X,  f22,  f23;  f42,  f43;  f63)  were 
found  to  satisfy  a boundary  value  differential  system  consisting  of  the  Clairaut  equations 


a2/"  + 6 aDf  + (60  - C)  f = R 

with  the  following  conditions  at  both  ends  of  the  interval  (0,  at): 

A0)  = f(0)  = 0 
Aftax)  + Bf,{ax)  = S(a,). 


(0) 


(10) 
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Hor*«  the  primes  denote  derivatives  with  respect  to  the  mean  radius  a,  and  A , B,  C arv 
constants  depending  on  the  order  of  approximation.  The  function 


W(*>)  = P0(ul/p0(u) 

is  the  ratio  between  the  density  py(a)  and  the  mean  density 


Po(°) 


p0u2<fo. 


til) 


(12) 


/f(u)  and  3(a)  are  known  functions  of  u that  depend  on  lower  order  approximations. 

lhe  specific  values  of  the  constants  and  the  form  of  the  functions  ap|x*anng  in 
equations  (9)  and  (10)  were  given  in  the  above  mentioned  papers  and  will  not  l>o  repro- 
duced here. 


EARTH  DENSITY  MODELS 

One  must  integrate  the  Clairaut  equation  to  obtain  the  ellipticity  of  the  equipotential 
surfaces  as  well  as  the  rotational  deformation  of  the  outermost  equilibrium  surface.  For 
this  purpose,  a density  profile  p^(r)  should  be  provided.  In  this  paper  we  plan  to  take  into 
consideration  recent  Earth  density  nuMlels,  to  integrate  the  t'lairaut  equation,  and  to  com- 
ixire  results  pertaining  to  the  exterior  potential  with  data  obtainable  from  satellite  geodesy. 

In  a recent  paper,  Lansano  and  Daley  131,  the  Hadden  and  Bullen  1969  HB1  Earth 
density  nuxlel  was  discussed  and  the  results  of  the  numerical  integration  were  compared 
with  earlier  work  by  James  and  Kopal  (4).  For  numerical  integration  purposes  we 
have  user!  hen*  the  1975  1066A  Earth  density  model  by  Gilbert  and  Dziewonski  and  the 
1977  QM2  nuxlel  by  Anderson  and  Hart. 

A brief  discussion  of  these  models  ap|M*ars  appropriate  and  will  lx*  provided  here  to 
enable  the  reader  to  understand  their  relative  merits  and  limitations  and  thereby  compre- 
hend the  validity  of  the  numerical  results. 

Bullen  has  provided  the  first  approximate  Earth  models  by  pioneering  work  in  which 
he  list'd  the  fact  that  the  travel  times  of  the  Ixxlily  waves  do  not  depend  upon  the  source 
mechanism  of  the  earthquake  but  are  functionals  of  structure  alone.  These  models  were 
achieved  by  relating  the  travel  time  with  the  travelled  angular  distance  and  by  making  use 
of  the  following  equations: 


«« 

Pq  3 


.p  , e. ; 


Pi)  di 


= »)P(Wo/** 


(13, 


3 


\ 


LANZANO  AND  DALKY 


where  k p,  a,  0 donate.  respectively,  the  incompressibility,  the  rigidity,  and  the  P and 
S seismic  velocities  at  depth  * from  the  surface  or  at  dutww  r^m  Uie  center  here  n 
is  a coefficient  depending  on  the  homogeneity  of  the  material  (see.  e.g..  Bull  (5|, 
pp.  227-240). 

These  early  Bullen  models  were  used  by  Kopal  [4]  to  determine  the  second  and 
fourth  harmonics  of  the  geopotential  from  hydrostatic  theory  alone  by  integrating  a 
second-order  approximation  of  the  Clairaut  equation. 

Improvements  to  these  early  models,  which  were  based  on  travel  time  data  alone, 
can  lie  obtained  by  making  use  of  the  observed  frequencies  of  the  Earth  vibrations  caused 
by  earthquakes.  These  vibrations  can  be  construed  as  being  the  superposition  of  the 
elastic-gravitational  normal  modes  of  the  Earth  that  are  excited  by  an  earthquake  Om 
can  construct  average  Earth  models  through  the  inversion  of  the  observed  eigenpcriods 
of  the  Earth.  Mathematically  the  problem  consists  of  finding  perturbations  to  a given 
Earth  model  so  that  the  differences  between  calculated  and  observed  eigenfrequenc.es 
can  be  minimized. 

The  normal  modes  of  oscillation  for  a nonrotating,  spherically  symmetric  Earth  an* 
of  two  kinds:  spheroidal  and  toroidal.  The  amplitudes  of  their  oscillations  are  represen- 
table, respectively,  as 


= ntMr)'T(0,0)  r + «>.*). 

n7*m  = (r)  r X VY7  (d ,<p). 


(14) 


These  oscillations  depend  on  three  parameters:  the  angular  order  £.  the  azimuthal  order 
m (with  -I  < m < I)  - which  are  related  to  the  degree?  and  order  m of  the  spherical 
harmonic  Y™  ((),</>)  - and  the  radial  number  n,  which  defines  the  overtone. 

The  azimuthal  urde,  m due,  hot  a, .fern  m the  differential  o,,„«ti<m»  I or 
ditions  so  that  for  each  land  n there  an*  21  + 1 normal  mode  eigenfunctions,  all  of  which 
belong  to  the  same  eigenfroquency;  together  they  form  n multiplot.  Dus  Ph^TJ'S 
called  degenernev  and  constitutes  the  main  difficulty  in  identifying  spectral  (leaks  o t 
eigen  frequencies  and  in  comparing  them  with  theoretical  calculations.  Perturbations,  smh 
annotation  ellipUcitv.  and  geographical  features  remove  the  degeneracy  of  a multiple*  and 
“a^bi  a sph  mu  tSet.  Ascertaining  the  splitting  parameter,  of  a multiset  would  be  a 
Zl  lead , ng  toth  e possibility  of  solving  for  singlets;  however,  the  theory  behind  such 
splitting  is  not  well  understood  and  consequently  has  not  been  formulated^  A«wdigl> . 
all  the  Earth  models  available  up  to  the  present  time  are  based  on  a nonrotating  ha  tl  . 

On  the  other  hand,  new  developments  in  instrumentation  allow  geophysicists  to  U wtr 
the  detection  threshold  for  magnitudes  in  the  low  frequency  modes  and  thus  mcreas. 
the  number  of  usable  seismic  events. 
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In  the  1969  HB1  Earth  model  by  Haddon  ami  Hullen  |6],  110  observations  of  |ieriods 
of  spheroidal  n.ul  toroidal  oscillations  were  taken  into  account  baaed  on  the  earthquakes 
of  May  1900  in  Chile  and  March  1904  in  Alaska.  Also,  for  the  first  time,  use  was  made 
of:  (U  the  revised  value  of  the  Earth’s  polar  moment  of  inertia:  lf,  0.3309  m,u'f , 
where  m,  5.970. 11127*!  is  the  Earth’s  total  mass  and  <i,  0371  km  is  the  Earth’s  mean 

radius  (i.e.,  the  radius  of  a sphere  of  equal  volume)  and  (2)  the  evidence  that  the  central 
density  pc  is  s 13  g/cma. 

In  the  preceding  as  well  as  in  many  other  free  oscillation  studies  of  the  Earth's  inter- 
ior, the  effect  of  absorption  u)H>n  the  eigen jieriods  of  the  Earth  has  been  ignored ; tins 
is  equivalent  to  assuming  that  the  Earth  is  perfectly  elastic.  The  actual  Earth  is  significantly 
anelestic  over  seismic  frequencies,  and  recent  research  by  Hart  et  al.  |7]  has  revealed 
that  a frequency-deiHMident  correction  of  the  order  of  1%  should  be  applied  to  the  normal 
mode  periods  to  eliminate  baseline  discre|vincies  between  body  wave  results  and  normal 
mode  results.  These  authors  adjusted  the  eigenperiods  of  their  older  C2  Earth  model  for 
attenuation  and  then  by  inversion  obtained  the  QM2  model  (1977).  Inclusion  of  the 
attenuation  term  tends  to  increase  the  values  of  seismic  velocities,  esjiecially  the  shear 
velocity. 

All  the  methods  mentioned  are  I vised  on  disjiersion  characteristics  of  traveling  waves, 
the  eigenfrequencies  of  normal  modes,  and  the  attenuations  of  both,  which  an*  essentially 
functionals  of  the  Earth  structure  alone.  A complementary  procedure  for  improving  upon 
the  mechanical  model  of  the  Earth  is  to  study  the  source  mechanism  of  an  earthquake  from 
the  knowledge  of  the  normal  mode  amplitudes-.  Using  both  methods  one  can  then  initiate 
an  iterative  process  of  successive  refinements  of  structure  and  mechanism. 

This  is  essentially  the  procedure  followed  by  Hilbert  and  Driewonski  |8]  in  obtaining 
their  1 066 A model  in  1975.  The  gross  Earth  data  were  compiled  from  1461  modes, 
whereas  the  source  mechanism  (or  moment  tensor)  was  retrieved  from  the  seismic  sjvvtra 
of  the  July  1970  Colombia  and  the  August  1963  Peru-Bolivia  earthquakes.  Tables  1.  2. 
and  3 represent  the  density  variations  with  depth  for  the  llBl,  1066A  and  QM2  models, 
resjiectively. 


NUMERICAL  INTEC, RATION  OE  THE  CLAIR AUT  EQUATIONS 
AND  RESULTING  EARTH  DEFORMATIONS 

We  have  used  central  difference  formulas  to  express  the  first  and  second  derivatives 
of  the  unknown  functions  at  the  various  pivotal  points  chosen  within  the  integration  range. 
Neglecting  thin!  and  higher  order  central  difference  terms  and  assuming  equal  intervals 
lietween  the  pivotal  points,  we  find  that  the  Clairaut  equation  (9)  takes  the  form  of  a 
three-term  set  of  difference  equations: 

*0.1  ♦ .4(i.2)f,  ♦ .4(i,3)/rlM  .4(1,4);  (i-  1.2 N).  (15) 

where  the  subscript  i refers  to  the  pivotal  point  at  which  the  functions  are  evaluated,  and 
the  .4 's  .lejiend  on  the  coefficients  of  the  original  equation  and  the  interval  sue. 
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Table  2 Values  of  Density  p,  in  the  Karth  Model  1066A 


11 

11 

37 

62 

88 

114 

139 

165 

190 

216 

242 

267 

293 

319 

344 

370 

395 

421 

421 

452 

484 

515 

546 

577 

609 

640 

671 

671 

706 

740 

775 

810 

844 

879 

913 

948 

983 

1,017 

1,052 

1,087 


3.372 

3.379 

3.387 

3.393 

3.402 

3.419 

3.443 

3.473 


3.600 

3.657 

3.712 

3.712 

3.764 

3.805 

3.850 

3.903 

3.961 

4.025 

4.106 

4.208 

4.208 

4.319 

4.403 

4.468 

4.511 

4.535 

4.537 

4.537 

4.542 

4.552 

4.567 

4.587 

4.610 


,121 
.156 
.191 
,225 
,260 
1,294 
1,329 
1,364 
1,398 
1,433 
1,468 
1,502 
1,537 
1,571 
1 .606 
1,641 
1.675 
1,710 
1.745 
1,779 
1,814 
1,849 
1,883 
1,918 
1,952 
1,987 
2,022 
2,056 
2.091 
2,126 
2,160 
2,195 
2,229 
2,264 
2,299 
2,333 
2,368 
2,403 
2,437 
2,472 
2,507 
2,541 


4.874 
4.892 
4.91 1 
4.930 
4.949 
4.969 
4.988 
5.008 
5.027 
5.046 
5.066 
5.086 
5.106 
5.127 
5.147 
5.169 
5.190 
5.212 
5.233 
5.255 
5.277 
5.298 
5.318 
5.338 
5.357 
5.374 
5.391 
5.406 


2,576 

2,610 

2,645 

2,680 

2,714 

2,749 

2,784 

2,818 

2,863 

2,887 

2,887 

2,958 

3,028 

3,099 


3.592 

3,663 

3,733 

3,803 

3,874 

3,944 

4,015 

4.085 

4,156 

4,226 

4,297 

4,367 

4,438 

4,508 

4,578 

4,649 

4,719 

4,790 

4,860 

4,931 

5,001 

5,072 


5.447 

5.460 

5.471 

5.483 

5.495 

5.506 

5.518 

5.528 

9.914 

10.028 

10.134 

10.235 

10.333 

10.427 

10.516 

10.603 

10.687 

10.772 

10.858 

10.946 

11.033 

11.116 

11.192 

11.265 

11.335 

11.406 

11.475 

11.541 

11.602 

11.660 

11.716 

11.769 

11.818 

11.863 

11.904 

11.944 

11.985 

12.026 

12.068 

12.110 


5,142 

12.153 

5,142 

13.021 

5,181 

13.031 

5,219 

13.045 

5,257 

13.060 

5,296 

13.074 

5,334 

13.088 

5,373 

13.102 

5,4 1 1 

13.116 

5,449 

13.131 

5.488 

13.145 

5.526 

13.160 

5,565 

13.175 

5,603 

13.190 

5,641 

13.205 

5.680 

13.220 

5,718 

13.235 

5,757 

13.251 

5.795 

13.267 

5,833 

13.284 

5,872 

13.300 

5,910 

13.316 

5,949 

13.331 

5,987 

13.345 

6,025 

13.357 

6,064 

13.368 

6,102 

13.377 

6,141 

13.385 

6,179 

13.393 

6,217 

13.399 

6,256 

13.406 

6,294 

13.411 

6,333 

13.418 

6,371 

13.421 
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Table  3 — Values  of  Density  p,  in  the  Earth  Model  QM2 


Depth  p I 

(km)  (g/cm3)  DePth  P Depth  p I Depth  p 
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Tables  4,  5,  and  6 refer  to  the  HB1,  1066A,  and  QM2  Earth  density  models,  respectively; 
they  furnish  the  deformation  coefficients  fQ , f2,  and  f6  for  the  equipotential  surfaces 
corresponding  to  intermediate  values  of  the  mean  radius  a.  The  number  N of  pivotal 
points  required  to  reach  the  mentioned  accuracy  was  851,  580,  and  607,  respectively. 

The  results  tabulated  at  the  chosen  values  of  the  depth  are  linear  interpolations  between 
the  pivotal  points. 


Table  4 — Deformation  Coefficients  for  the 
Equipotential  Surfaces,  Earth  Model  HB1 


Depth 

(km) 

-106  • f0 

~102  • f 2 

+ 10&  . 

00 

o 

1 

0 

0.99421 

0.22298 

0.44766 

0.96569 

15 

0.99147 

0.22268 

0.44630 

0.96106 

300 

0.94109 

0.21695 

0.42139 

0.87871 

350 

0.93259 

0.21596 

0.41724 

0.86551 

650 

0.88288 

0.21013 

0.39312 

0.79141 

900 

0.84214 

0.20522 

0.37324 

0.73257 

1200 

0.79332 

0.19918 

0.34890 

0.66167 

1500 

0.74495 

0.19302 

0.32400 

0.58968 

1800 

0.69804 

0.18684 

0.29891 

0.51695 

2100 

0.65431 

0.18089 

0.27448 

0.44533 

2400 

0.61654 

0.17559 

0.25242 

0.37940 

2878 

0.58072 

0.17041 

0.23085 

0.31418 

3000 

0.57747 

0.16994 

0.22905 

0.30954 

3300 

0.57030 

0.16888 

0.22508 

0.29935 

3600 

0.56416 

0.16797 

0.22167 

0.29071 

3900 

0.55893 

0.16719 

0.21876 

0.28338 

4200 

0.55451 

0.16652 

0.21628 

0.27720 

4500 

0.55084 

0.16597 

0.21422 

0.27205 

4800 

0.54786 

0.16552 

0.21253 

0.26785 

5121 

0.54544 

0.16516 

0.21114 

0.26440 

5400 

0.54385 

0.16492 

0.21023 

0.26215 

5700 

0.54260 

0.16473 

0.20952 

0.26038 

6000 

0.54182 

0.16461 

0.20907 

0.25927 

6300 

0.54149 

0.16456 

0.20888 

0.25880 

6371 

0.54147 

0.16455 

0.20887 

0.25878 
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Table  7 shows  the  surface  values  (i.e.,  those  values  for  a = a1)  of  the  f's  and  of  the 
coefficients  K appearing  in  the  spherical  harmonic  expansion  of  the  potential  according 
to  the  formula 


V0(r,  cos0,  a1)  = ■ 


Gm  i 


1 + 


^ PytootB) 

j- 1 


(16) 


Also  shown  is  the  value  of  the  ellipticity, 

r(ax,  0)  - r(ax,  1) 


e = 


r(a i,  0) 


(17) 


and  of  the  ratio  (C  -A  )/C,  where  A is  the  moment  of  inertia  with  respect  to  any  barycentric 
axis  in  the  equatorial  plane  and  C the  moment  of  inertia  with  respect  to  the  polar  axis. 

This  table  is  a comparison  of  four  Earth  density  models  and  includes  the  results  for  the 
Bullen  model,  which  required  only  387  pivotal  points  for  reaching  the  same  accuracy. 


Table  7 — Comparison  of  Earth  Density  Models* 


Parameter 

QM2 

1066A 

CQ 

X 



Bullen  (1940) 

fo  (al) 

-0.99239  • 10-6 

-0.99376  • lO-6 

-0.99421  ■ 10-6 

-0.10098  • 10-5 

f 2 (al) 

-0.22278  • 10-2 

-0.22293  • lO-2 

-0.22298  • 10-2 

-0.22473  • 10-2 

f 4 (fll) 

+0.44648  • lO"6 

+0.44762  • 10-6 

+0.44766  • lO'5 

+0.45202  • 10-6 

f 6 (°l) 

-0.96148  • 10-6 

-0.96520  ' lO"8 

-0.96569  • lO"8 

-0.97129  ■ lO"8 

^2  (°1 ) 

-0.10735  • 10-2 

-0.10751  • lO'2 

-0.10756  • lO"2 

-0.10929  • lO'2 

K4  («x) 

+0.29625  • 10"5 

+0.29763  • lO"6 

+0.29776  • lO"6 

+0.30495  ' lO-6 

^6  (al) 

-0.11199  • 10-7 

-0.11284  • 10-7 

-0.11293  • lO"7 

-0.11598  • lO’7 

e-i 

299.8 

299.6 

299.6 

297.2 

A = Ie 

0.80087  1 036 

0.80168  • 1036 

0.80212  • 1036 

0.80893  • 1036 

C = Ip 

0.80347  • 1036 

0.80429  • 1036 

0.80473  • 1036 

0.81159  • 1036 

(C-A)IC 

0.3239  • 10-2 

0.3242  • 10-2 

0.3243  • lO’2 

0.3274  • 1C’2 
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Wo  romark  here  that  tho  geodetic  valuos  obtained  from  satollito  moasuromonts  lie 
between  the  data  for  the  Bullen  and  the  HB1  models.  Such  deviations  can  lx-  assumed  to 
represent  a measure  of  the  role  the  elastic  energy  plays  in  the  rotational  deformation  of 
the  Karth. 


Notice  also  the  apparent  tendency  of  the  three  new  models  to  converge  toward  limi- 
ting values  of  the  parameters,  a trend  that  augurs  well  for  the  future  knowledge  of  a real 
Karth  model. 


DEFORMATIONS  OF  JUPITER  AND  SATURN 

Jupiter  ami  Saturn  comprise  about  90%  of  the  total  mass  of  the  planetary  system; 
therefore,  the  great  interest  which  in  recent  years  has  arisen  toward  developing  a density 
model  for  their  interior  and  ascertaining  the  oblateness  of  their  exterior  shape  is  quite 
understandable. 

Pioneering  work  on  Jupiter’s  internal  composition  goes  back  to  Wildt  [10)  and 
Ramsey  1 1 1 1 . Their  theory  was  further  elaborated  by  DeMarcus  |12|.  Opik  |13|, 
and  Peebles  1 1 4 1 . More  recently,  Podolak  and  Cameron  1151,  Zharkov  et  al.  1161, 
and  Slattery  |17|  have  contributed  to  the  topic. 


In  the  above  mentioned  works,  the  Jupiter  interior  is  conceived  as  a hydrogen  and 
helium  envelope  approximately  in  the  solar  mixture  surrounding  a core  of  heavy  elements 
By  solar  mixture  we  mean  that  the  hydrogen  mass  fraction  is  0.78  and  the  helium  mass 
fraction  is  0.22.  Within  this  gaseous  envelope,  one  can  distinguish  primarily  three  layers 
beginning  with  a very  low  density  region  that  can  be  treated  adequately  with  perfect  gas 
laws,  a transition  region  of  gas  mixture  in  a molecular  state  where  the  Van  der  Waals 
equation  can  be  applied,  anil  a metallic  hydrogen  region  for  which  the  equation  of  state 
is  also  known,  lhc  various  models  which  have  been  considered  differ  primarily  because 
of  the  equation  of  state  adopted  for  the  transition  region. 

Saturn’s  models  have  evolved  concomitantly  with  Jupiter’s.  However,  since  Saturn 
is  smaller  than  Jupiter,  more  of  the  planet  ’s  interior  is  expected  to  be  in  the  molecular 
state,  and  thus  more  uncertain  in  its  composition. 

DeMarcos  [121  made  use  of  the  second-order  theory  of  hydrostatic  equilibrium 
as  developed  by  DeSitter  1 181  and  recast  it  into  a form  which  is  very  suitable  for 
ascertaining  the  density  profile  of  a planet  when  only  the  equation  of  state  is  known. 

Notwithstanding  claims  made  later  by  Slattery  |17),  Zharkov  (16)  used  only 
a second-order  theory  in  developing  Jupiter’s  density  profile.  To  l»e  more  specific,  in 
1968  Zharkov  1 19)  did  develop  a third-order  equilibrium  theory  which  appeared  in  English 
translation  in  Soviet  Astronomy  (1970);  however,  he  limited  himself  to  obtaining  integral 
relationships  basically  equivalent  to  those  already  achieved  by  lanzano  in  1962  ( 201  • The 
complicated  form  assumed  by  Zharkov’s  equations  makes  it  unlikely  that  he  could  have 
used  such  theory  to  obtain  Jupiter’s  density  profile. 
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Slattery  (17)  developed  two  Jupiter  models  in  accordance  with  which  core  density 
reached  values  of  13  and  22  g/cma.  In  developing  these  models,  he  used  a fourth-order 
theory  elaborated  by  Hubbard  et  al.  121).  This  theory  is  objectionable  not  only  from 
a theoretical  point  of  view  because  no  use  is  made  of  the  equipotential  surfaces,  but  also 
from  a computational  viewpoint  because  of  the  doubtful  convergence  of  the  expressions 
therein  considered. 

Because  of  these  drawbacks,  in  the  present  work  we  use  the  well-tested  lJeMarcus 
density  profile  (Table  8)  in  conjunction  with  the  third-order  theory  of  Kopa!  and  Iamzano 
to  obtain  more  accurate  values  for  the  deformation  coefficients  of  the  equipotential  sur- 
faces and  for  the  coefficients  of  the  harmonics  in  the  exterior  potential.  This  third -order 
theory  is  the  analytic  elaboration  of  our  already  mentioned  1962  results,  whereby  we  were 
successful  in  transforming,  by  means  of  differentiation  and  elimination  operations,  a rather 
intricate  system  of  integral  equations  into  a simple  system  of  second-order  ordinary  dif- 
ferential equations  with  boundary  conditions. 

The  initial  data  we  used  for  Jupiter  and  Saturn  have  been  taken  from  Brouwer  anil 
Clemence  (22)  and  art'  summarized  in  Table  9. 
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Table  9 Data  of  Ihrouwer  anil  (.'lenience  for  Major  Planets* 


Parameter 

Jupiter 

Saturn 

m, 

1.902  1030g 

0.5694  • 103°g 

al 

6.9861  • 10®  cm 

5.763  • 10® cm 

Po  (°t> 

1.33  g/cm3 

0.71  g/cm3 

Po  (9) 

30.84  g/cm3 

15.62  g/cm3 

T 

9.87  hr 

10.41  hr 

1.768317  • 10"*  sec*1 

1.676589  10*  4 six'-1 

q 

0.28004  • l0-> 

0.47207  10*1 

•Fumlanu'nul  pammolors  uw'il  in  the  integration. 


Results  of  the  Pioneer  10  ami  1 1 missions  nave  rise  to  revised  data  for  Jupiter  (see 
Anderson  et  al„  |23|),  primarily  for  its  radius  o,  7.14  • 10 9 em,  and  its  rotational  per- 
iod T 0.925  hours,  whereby  q 0.0200.  We  have,  however,  used  the  older  sot  of  data 
not  only  for  the  sake  of  com|mrison  with  Kopal’s  results,  hut  also  heeause  we  do  not  be- 
lieve that  the  newer  data  wmdd  have  changed  the  results  appreciably. 

Our  results  are  summarized  in  Table  10.  To  achieve  accurate  solutions,  we  had  to  take 
284  pivotal  points  in  the  integration  of  Jupiter’s  deformations,  and  401  points  in  the  case 
of  Saturn’s.  The  ellipticity,  « , is  0.004433  for  Jupiter  and  0.097847  for  Saturn. 

We  compare  our  results  with  the  following: 

• Observational  results  of  Brouwer  and  Clemence  and  of  Anderson  et  al.  (Table  11); 

• 1963  results  by  James  and  Ko|»al  (Table  12),  where  the  DeMarcus  profile  and  a 
second -order  theory  were  used;  and 

• results  in  1977  by  Slattery  (Table  13);  he  used  a fourth-order  approximation,  the 
new  data,  and  a different  density  profile. 
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THE  LINEARIZED  NAVIER-STOKES  EQUATION  FOR  AN  ELASTIC  EARTH 

To  go  one  step  further  and  consider  a realistic  model  of  the  Earth,  we  must  take  into 
account  the  effect  of  the  elastic  forces.  For  this  purpose,  we  assume  that  the  additional 
stress  field  Ty  and  the  corresponding  displacement  field  n|(  measured  from  the  reference 
state,  art'  related  according  to  tin*  relationship 
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Table  10  — Surface  Values  for  Major  Planets 


| Jupiter 

Saturn 

Parameter 

fo  <*»i ) 

-0.39517  • 10-3 

-0.94731  • 10*3 

f 2 (®l) 

-0.44701  • 10'1 

-0.69585  • 10- 1 

f 4 <«1> 

+0.20037  • lO-2 

+0.57156  • lO'2 

f 6 (°l) 

-0.91101  • 10-4 

-0.42411  • 10-3 

K2  (aj) 

-0.15013  • 10- 1 

-0.18465  • lO'1 

*4  («»l) 

+0.68165  • lO"3 

+0.15885  • lO'2 

^6  (al ) 

-0.39233  • 10-4 

-0.15591  • 10-3 

e-l 

15.52 

10.22 

A =Ie 

0.23691  • 1040 

0.39762  • 1030 

C = 'p 

0.25068  • 1040 

0.43157  • 1039 

(C-A)IC 

0.5494  • 10"1 

0.7867  • 10- 1 

Table  11  — Results  of  Brouwer  and  Clemence  and  of  Anderson 


Parameter 

Jupiter 

Saturn 

*2 

*4 

e 

* 

-0.147066  • lO’1 

0.674666  • 10-3 

rr-~  - 0.0651890 
15. o4 

** 

-0.1472  • 10- 1 

0.65  lO"3 

* 

-0.166733  • 10"1 

0.102933  • lO*2 

10,21  ' 0 0979432 

•Brouwer  and  Clemence  (1961) 
••Anderaon  et  al.  (1974) 

|1 


Table  13-Data  of  Slattery  (1977) 


Fourth-order  approximation,  new  data,  and  a different 
density  profile  were  used. 


This  formula,  in  which  use  has  been  made  of  the  summation  convention,  is  valid  for  a 
perfectly  elastic  and  isotropic  medium.  Here  X and  p are  the  Lame  elastic  parameters; 

H is  also  known  as  the  rigidity.  These  two  parameters  are  related  to  the  incompressibility 
or  bulk  modulus  by  the  relationship 


k = \+jn.  (19) 
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Within  the  validity  of  the  infinitesimal  deformation  theory  we  can  safely  assume  that 
equation  (18)  applies  to  the  coordinates  that  the  points  of  the  medium  had  before  the 
occurrence  of  deformation.  This  assumption  cannot,  however,  be  made  for  the  initial 
stress  field;  as  a consequence,  the  total  stress  field  at  the  undeformed  points  must  be 
written  as 


°ii  = Tit  uk  + r0  . (20) 

With  respect  to  an  inertial  frame,  the  equations  representing  the  deviation  from  the  ref- 
erence state  are 


p — -pit  +l3t. 

dt3  9*,  dxj 


(21) 


p is  the  mass  density  and  is  the  sum  of  two  terms: 

P = Po  + Pt  • 

where  p,,  the  change  due  to  the  displacement  field  u,  is  given  by  the  continuity  equation 


*'  S' 


(22) 


V = v0  + V,  is  the  total  potential,  where  V,  satisfies  the  Poisson  equation 

V*vl  = -4»Gpj  . (23) 

The  components  of  gravity  will  then  appear  to  l>e 

= *0,i  + *i,i  . (2-1) 


where  g,  , = 


Neglecting  the  product  Pjffj  ,,  since  it  is  of  the  second  order  in  the  displacements, 
and  considering  that 

3Ty  _ 3/>o_ 

djr*  " A**  5</  • 

we  arrive  at  the  following : 

da«,  ft  0r.. 

P~dt*  = Po*'’'  + Pl*°’'  + + ' (2M 
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Note  that  equations  (1),  (20)  and  (21)  have  been  used  in  obtaining  it.  The  last  term  in 
equation  (25),  when  expanded,  yields 


3r0  _ 

/3X 

duk 

+ x \ 

6 tut*2"* 

+ 92-A 

/dp,. 

. 

dXj 

\dxj 

dxk 

dxjdxj 

v \ dxjdxj 

3x,3x;7 

1 + 

3 Xj 

u 

3 J 

The  reduction  of  the  previous  expressions  to  vectorial  notation  can  be  accomplished  by 
taking  into  account  the  following  facts: 


duk 

1)  t — = V • u is  the  dilatation  A; 
°xk 


2) 


V(V-  u); 


d2Uj 

fa. qx.  is  the  lth  component  of  the  gradient  of  the  dilatation, 


l 

t 


‘ 


32u,- 

3)  fa. fa.  *s  ^e  component  of  the  vector  Laplacian  V2u,  which  is  known  to  be 
representable  according  to  the  vector  identity  V2u  = V(  V-  u)  - V X V Xu; 


4)  Finally  one  can  prove  (see  the  Appendix  for  details)  that  the  expression 

3p  / duj  dUj\ 
dxj  (ax,  + dxj 

is  the  ith  component  of  the  vector 

(VP)(  V • u)  + v(u  • Vp)  + V X (u  X Vp)  -uV2p.  (27) 

Replacing  the  above  vector  quantities  within  equations  (25)  and  (26),  one  gets 

d2  u 

P df2"  = P°9l  + Pl9o  + V(Pou  • 9o>  + (*  + 2p)  V(V  • u)  (28) 

-HV  X V X u 4-  (VX  + Vp)(V  ■ u)  + V(u  ■ Vp)  + V X (u  X Vp)  - uV2p  . 


This  is  essentially  the  Navier-Stokes  equation  when  allowance  is  made  for  the  variation 
of  the  material’s  elastic  parameters  and  when  quadratic  terms  in  the  displacements  are 
dropped.  X,  p,  and  p0  are  supposed  to  be  known  functions  of  the  radial  distance.  As  a 
matter  of  fact,  it  easily  follows  from  equations  (13)  and  (19)  that  the  Lame  parameters 
are  related  to  the  P and  S seismic  velocities  according  to 
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p0  a2  = X + 2p;  p0/32  = p. 

Consequently,  their  expressions  as  functions  of  the  radial  distance  or  depth  can  be  obtained 
from  the  density  and  the  velocity  of  seismic  waves.  Their  values  are  provided  by  the  Earth 
model  tables. 

g0  and  g1  are  obtainable  from  V0  and  Vj,  respectively,  via  the  Poisson  equations  (4) 
and  (23);  on  the  other  hand,  px  is  obtainable  from  the  continuity  equation  (22). 

The  acceleration  of  u with  respect  to  an  inertial  frame  is  given  by 


d2u  d2u  du  dco 

df2  d(2  dt  dt 


X u + (gj  • u)  CO  - CJ2U, 


(29) 


where  the  symbol  of  partial  derivative  refers  to  the  variation  of  a vector  with  respect  to 
the  rotating  frame.  If  we  want  to  ascertain  a permanent  deformation,  the  vector  u, 
measured  with  respect  to  a rotating  Earth,  should  be  independent  of  time;  this  means 
that 


d2u  _ du  _ 
dt2  ~ dt 

Also,  in  the  case  of  a constant  angular  rate  and  of  no  variation  in  the  position  of  the  ro- 
tational axis,  one  should  have 


d co 
dt 


= 0. 


Therefore,  we  are  left  with  the  expression 

d2u 

= (co  • u)  to  - co2u.  (30) 

dt2 

Using  equation  (30),  we  find  that  the  left-hand  side  of  equation  (28)  becomes 

p0  [(co  • u)  co  - to2u]  , (31) 

where  the  term  ptu  has  been  neglected  because  it  is  of  the  second  order  in  the  displace- 
ments. 


When  the  left-hand  side  is  replaced  by  equation  (31),  equation  (28)  represents  the 
fundamental  relation  which  yields  the  perturbation  displacement  u.  One  should  seek 
solutions  to  this  equation  in  the  form  of  spheroidal  and  toroidal  deformations  as  expres- 
sed by  equation  (14). 
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DISCUSSION  OF  THE  RESULTS 

We  have  developed  a general  equilibrium  theory  which  has  two  functions: 

(1)  It  yields  the  fundamental  state  of  deformation  for  a rotating  body  in  the  presence 
of  hydrostatic  forces;  this  can  be  accomplished  by  solving  the  Clairaut  equation  with  boun- 
dary conditions;  and 

(2)  It  accounts  for  the  additional  deformations  attributable  to  elastic  forces  as  per- 
turbations to  the  previously  obtained  fundamental  mode;  the  perturbation  equation  in 
question  is  a linearized  version  of  the  Navier-Stokes  equation. 

For  the  Earth  we  have  numerically  solved  the  Clairaut  equation  using  various  density 
profiles  and  have  compared  our  results  with  well  established  satellite-obtained  data  per- 
taining to  the  second,  fourth,  and  sixth  harmonics  of  the  geopotential.  The  observed 
discrepancies  are  primarily  due  to  the  presence  of  shear  stresses  and  convection  currents 
in  the  Earth’s  mantle.  To  remedy  this  situation,  one  must  solve  the  Navier-Stokes  per- 
turbation equation;  this  we  plan  to  do  in  the  near  future  as  our  next  task. 

In  the  case  of  the  two  giant  planets,  our  results  represent  an  improvement  on  Kopal’s 
previous  work,  which  was  based  on  the  second-order  Clairaut  equation  and  the  same 
DeMarcus  density  model.  Our  results  compare  favorably  with  recent  work  by  Hubbard 
and  Slattery,  who  use  a different  density  distribution,  primarily  because  of  our  more 
advanced  mathematical  formulation:  preliminary  analytical  work  has  allowed  us  to 
simplify  the  original  equilibrium  conditions  which  as  used  by  those  authors  were  expressed 
by  a complicated  integro-differential  relationship;  we  obtained  a well-posed  boundary  value 
problem  consisting  of  the  Clairaut  equation.  This  feature  of  simplicity  is  responsible  for 
obtaining  more  accurate  solutions  more  rapidly. 

A preliminary  version  of  these  results  was  presented  orally  at  the  Spring  Meeting  of 
the  American  Geophysical  Union  in  Miami  Beach,  April  1978. 
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Appendix 


MATHEMATICAL  PROOF  OF  A VECTOR  IDENTITY 
From  the  definitions 


(u  X Vp)„ 


~ fnjk 


(V  X u)j  = 


du^ 

t ,m"  **m' 


ami 


where  e is  the  completely  antisymmetric  third-rank  tensor  (also  known  as  the  permutation 
symbol  or  Levi-Civita  density),  one  gets  the  following: 

IV  X („  X V(i>),  = ^ 

- ( ^u)  dp  . d2p 

dx*  U>  dx*dx 

— ^Ui  dp  ^uj  dp  d2p  d2p 

dxk  dxk  dxj  dx,  U|  dxfcdxfc  Uj  dx,dx;  * ^ 

Here  use  has  been  made  of  the  property  that  f is  antisymmetric  for  an  interchange  of  any 
pair  of  indices,  whereby 


enjk  fjkn  • 

also  the  tensor  identity, 

timnfjkn  ~ ^ij^mk  ~ ^ik^mj' 
was  used  — see,  e.g.,  Harris  (24|  pp.  10-13. 

Using  equation  (1  A),  one  can  express  the  ith  component  of  the  vector  represented  by 
equation  (27)  as 
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